Chapter 3 Week 3: Thematic maps in R

## Warning in Sys.setlocale("LC_ALL", "English"): OS reports request to set
## locale to "English" cannot be honored
## [1] ""

3.1 Intro and recap

Last week we showed you fairly quickly how to create maps of spatial point patterns using leaflet and we also introduced the tmap package for thematic maps. Besides doing that we introduced a set of key concepts we hope you have continued studying over the week. We also discussed the sf package for storing spatial objects in R.

This week we will carry on where we left the session last week. In the presentations last week we introduced various kind of thematic maps and in our lecture this week we discuss in detail issues with choropleth maps. So the focus of today’s lab is going to be around thematic maps and some of the choices we discussed in our presentation last week and also this week.

We will also introduce faceting and small multiples, which is a format for comparing the geographical distribution of different social phenomena. For this session we will be using the spatial object that you created last week and complement it with additional information from the census. So first of all you will have to rerun the code you used to create the “manchester_lsoa” sf object. Apart from doing so, you want to start your session loading the libraries you know for sure you will need:

You may not remember all of what you did to generate that file so let’s not waste time and just cut and paste from below (but try to remember what each of the lines of code is doing and if you are not clear look at the notes from last week). Imagine you had to do all of this again by pointing and clicking in a graphical user interface rather than just sending the code to the console! As you will see time and time again, code in the end is a much more efficient way of talking to a computer.

## Reading layer `england_lsoa_2011' from data source `/Users/reka/Dropbox (The University of Manchester)/crimemapping_textbook_bookdown/data/BoundaryData/england_lsoa_2011.shp' using driver `ESRI Shapefile'
## Simple feature collection with 282 features and 3 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 378833.2 ymin: 382620.6 xmax: 390350.2 ymax: 405357.1
## epsg (SRID):    NA
## proj4string:    +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs
## Warning: Column `code`/`LSOA.code` joining factors with different levels,
## coercing to character vector

You may not want to have to go through this process all the time. One thing you could do is to save the “manchester_lsoa object” as a physical file in your machine. You can use the st_write() function from the sf package to do this. If we want to write into a shapefile format we would do as shown below:

Notice how four files have appeared in your working directory, in your “BoundaryData”" subdirectory or whatever you called it. Remember what we said last week about shapefiles, there are a collection of files that need to be kept together.

If you wanted to bring this shapefile back into R at any future point, you would only need to use the st_read() function.

Before we carry on, can you tell What is different between manchester_lsoa.shp and manchester_crime_lsoa.shp? Think about it.

3.2 Creating choropleth maps

The tmap package was developed to easily produce thematic maps. It is inspired by the ggplot2 packaged and the layered grammar of graphics. It was written by Martjin Tennekes a Dutch data scientist. There are a number of vignettes in the CRAN repository and the GitHub repo for this package that you can explore. GitHub is a collaborative website used by software developers and data scientist, also contains a useful readme section with additional resources to familiarise yourself with this package. Each map can be plotted as a static map (plot mode) and shown interactively (view mode) as we briefly saw last week. We will start by focusing on static maps.

Every time you use this package you will need a line of code that specifies the spatial object you will be using. Although originally developed to handle sp objects only, it now also has support for sf objects. For specifying the spatial object we use the tm_shape() function and inside we specify the name of the spatial object we are using. On its own, this will do nothing apparent. No map will be created. We need to add additional functions to specify what we are doing with that spatial object. If you try to run this line on its own, you’ll get an error telling you you must “Specify at least one layer after each tm_shape”.

The main plotting method consists of elements that we can add. The first element is the tm_shape() function specifying the spatial object, and then we can add a series of elements specifying layers in the visualisation. They can include polygons, symbols, polylines, raster, and text labels as base layers. We will add a polygon using tm_polygon(). As noted, with tmap you can produce both static and interactive maps. The interactive maps rely on leaflet. You can control whether the map is static or interactive with the tmap_mode() function. If you want a static map you pass plot as an argument, if you want an interactive map you pass view as an argument. Let’s create a static map first.

## tmap mode set to plotting

Given that we are not passing any additional arguments all we are getting is a map with the shape of the geographies that we are representing, the census LSOAs for Manchester city. We can, however, ask R to produce a choropleth map by mapping the values of a variable in our data table using colour. In tmap we need to denote our variables between quotes. The first argument we pass then would be the name of the variable we want to visualise. If you remember we have a count for crimes (“count”), so let’s visualise that by creating a thematic map.

Notice how this map is different from last week. What do you think the reason for this is? You may remember last week we used one additional argument style specifying the classification method we were going to use. If you remember we used quantiles. We will in a second look at how a map of the counts of crime looks different when we use different classification systems. But before we get to that, let’s think about aesthetics a bit more.

We have been using tm_polygons() but we can also add the elements of a polygon map using different functions that break down what we represent here. In the map above you see the polygons have a dual representation, the borders are represented by lines and the colour is mapped to the intensity of the quantitative variable we are representing. With darker colours representing more of the variable, the areas with more crimes. Instead of using tm_polygon() we can use the related functions tm_fill(), for the colour inside the polygons, and tm_borders(), for the aesthetics representing the border of the polygons. Say we find the borders distracting and we want to set them to be transparent. In that case we could just use tm_fill().

As you can see here, the look is a bit cleaner. We don’t need to get rid of the borders completely. Perhaps we want to make them a bit more translucent. We could do that by adding the border element but making the drawing of the borders less pronounced.

The alpha parameter that we are inserting within tm_borders() controls the transparency of the borders, we can go from 0 (totally transparent) to 1 (not transparent). You can play around with this value and see the results.

Notice in the last few maps we did not have to specify whether we wanted the map to be static or interactive. When you use tmap, R will remember the mode you want to use. So once you specify tmap_mode(“plot”), all the subsequent maps will be static. It is only when you want to change this behaviour that you would need another tmap_mode call.

Notice as well that the legend in this map is (a) not very informative and (b) located in a place that is less than optimal, since it covers part of the map. We can add a title within the tm_fill to clarify what count is and we can use the tm_layout() function to control the appearance of the legend. This later function tm_layout allows you to think about many of the more general cosmetics of the map.

## tmap style set to "col_blind"
## other available styles are: "white", "gray", "natural", "cobalt", "albatross", "beaver", "bw", "classic", "watercolor"

3.3 Producing small multiples to compare the effect of different classification systems

For comparing the effects of using different methods we can use small multiples. Small multiples is simply a way of reproducing side by sides similar maps for comparative purposes. To be more precise small multiples are sets of charts of the same type, with the same scale, presented together at a small size and with minimal detail, usually in a grid of some kind. The term was at least popularized by Edward Tufte, appearing first in his Visual Display of Quantitative Information in 1983.

There are different ways of creating small multiples with tmap as you could see in the vignettes for the package, some of which are quicker but a bit more restricted. Here we are going to use tmap_arrange(). With tmap_arrange() first we need to create the maps we want and then we arrange them together.

Let’s make four maps, each one using a different classification method: Equal interval, Natural breaks (Jenks), Quantile, and Unclassed (no classification).

For each map, instead of visualising them one by one, just assign them to a new object. Let’s call them map1, map2, map3 and map4.

So let’s make map1. This will create a choropleth map using equal intervals as discussed earlier in the presentation:

Now create map2, with the jenks method often preferred by geographers:

Now create map3, with the quantile method often preferred by epidemiologists:

And finally make map4, an unclassed choropleth map, which maps the values of our variable to a smooth gradient.

Notice that we are not plotting the maps, we are storing them into R objects (m1 to m4). This way they are saved, and you can call them later, which is what we need in order to plot them together using the tmap_arrange() function.

So if you wanted to map just map3 for example, all you need to do, is call the map3 object. Like so:

But now we will plot all 4 maps together, arranged using the tmap_arrange() function. Like so:

## Some legend labels were too wide. These labels have been resized to 0.42, 0.42, 0.42, 0.42. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

3.3.1 Homework 1

Discuss which of these classification methods gives you the best visualisation of crime in Manchester city. Insert the produced maps in your answer.**

3.4 Using graduated symbols

The literature on thematic cartography highlights how counts are best represented using graduated symbols rather than choropleth maps. So let’s try to go for a more appropriate representation. In tmap you can use tm_symbols for this. We will use tm_borders to provide some context.

First thing you see is that we loose the context (provided by the polygon borders) that we had earlier. The border.lwd argument set to NA in the tm_symbols() is asking R not to draw a border to the circles. Whereas tm_borders() brings back a layer with the borders of the polygons representing the different LSOAs in Manchester city. Notice how I am modifying the transparency of the borders with the alpha parameter.

3.5 Bringing additional census data in

Last week you learned how to obtain crime data from the police UK website and you also developed the skills to obtain shapefiles with the boundaries for the UK census geographies. Specifically you learnt how to obtain LSOAs boundaries. Then we taught you how to join these data tables using dplyr. If you open your “manchester_lsoa” object you will see that at the moment you only have one field in this dataframe providing you with statistical information. However, there is a great deal of additional information that you could add to these data frame. Given that you are using census geographies you could add to it all kind of socio demographic variables available from the census.

You may want to watch this 4 minute video to get a sense for how to obtain the data. If you don’t have headphones make sure you read this brief tutorial before carrying on. We are going to get some data for Manchester city LSOAs. Let me warn you though, the census data portal is one of the closest things to hell you are going to come across on the internet. Using it will be a good reminder of why point and click interfaces can suck the life out of you.

From the main Infuse portal select the 2011 census data then when queried pick selection by geography:

Expand the local authorities and select Manchester. Expand Manchester and select LSOAs:

At the bottom of the page click in Add and then where it says Next. Now big tip. Do not press back in your browser. If you need to navigate back once you get to that point use the previous button at the bottom of the screen. You will regret it if you don’t do this.

Now you will need to practice navigating the Infuse system to generate a data table that has a number of relevant fields we are going to use today and at a later point this semester. I want you to create a file with information about: the resident population, the workday population, and the number of deprivation households. This will involve some trial and error but you should end up with a selection like the one below:

Once you have those fields click next to get the data and download the file. Unzip them and open the .csv file in Excel. If you view the data in Excel you will notice it is a bit untidy. The first row has no data but the labels for the variable names and the second row has the data for Manchester as a whole. We don’t need those rows. Because this data is a bit untidy we are going to use read_csv() function from the readr package rather than the base read.csv function.

## Warning: Missing column names filled in: 'X14' [14]
## Parsed with column specification:
## cols(
##   CDU_ID = col_double(),
##   GEO_CODE = col_character(),
##   GEO_LABEL = col_character(),
##   GEO_TYPE = col_character(),
##   GEO_TYP2 = col_character(),
##   F996 = col_character(),
##   F997 = col_character(),
##   F998 = col_character(),
##   F999 = col_character(),
##   F1000 = col_character(),
##   F1001 = col_character(),
##   F2384 = col_character(),
##   F323339 = col_character(),
##   X14 = col_logical()
## )

Notice that even all the variables that begin with “f” are numbers they have been read into R as characters. This is to do with the fact the first two lines do not represent cases and do have characters. R is coercing everything into character vectors. Let’s clean this a bit.

First we will get rid of the first two rows. In particular we will use the slice() function from dplyr. We can use slice to select cases based on row number. We don’t need the first two rows so we can select rows 3 to 284.

There are also fields that we don’t need. We only need the variables beginning with F for those have the information about population and deprivation, and the “GEO_CODE” tag which will allow us to link this table to the “manchester_lsoa” file.

We also want to convert the character variables into numeric ones, whilst preserving the id as a character variable. For this we will use the lapply function. This is a convenient function that will administer a function to the elements we pass as an argument. In this case we are asking to apply the as.numeric() function to the columns 2 to 9 of the “census_lsoa_m”" data frame. This is turning into numeric all those character columns.

The only problem we have now is that the variable names are not very informative. If you look at the metadata file that came along you can see that there is a key there to understand what these variables mean. We could use that information to create more meaningful names for the variables we have. We will use the rename() function from the dplyr package to do the renaming:

The rename function takes as the first argument the name of the dataframe. Then for each variable you want to change you write down the new name followed by the old name. Now that we have the file ready we can link it to our “manchester_lsoa” file using code we learnt last week. We use again the left_join() function to add to the “manchester_lsoa” dataframe the variables that are present in the “census_lsoa_m”. The first argument in the function is the name of the dataframe to which we want to add fields, the second argument the name of the dataframe from which those fields come, and then you need to specify using “by” the name of the variables on each of these two dataframes that have the id variable that will allow us to ensure that we are linking the information across the same observations.

And there you go… Now you have a datafile with quite a few pieces of additional information about LSOAs in Manchester. The next step is to use this information.

3.6 Computing and mapping crime rates

Ok, so now we have a field that provides us with the number of crimes and two alternative counts of population for each LSOA in Manchester in the same dataframe. We could compute the rate of crime in each using the population counts as our denominator. Let’s see how the maps may compare using these different denominators.

But first we need to create new variables. For this we can use the mutate() function from the dplyr package. This is a very helpful function to create new variables in a dataframe based on transformations or mathematical operations performed in other variables within the dataframe. In this function, the first argument is the name of the data frame, and then we can pass as arguments all new variables we want to create as well as the instructions as to how we are creating those variables.

First we want to create a rate using the usual residents, since crime rates are often expressed by 100,000 inhabitants we will multiply the division of the number of crimes by the number of usual residents by 100,000. We will then create another variable, “crimr2”, using the workday population as the denominator. We will store this new variables in our existing “manchester_lsoa” dataset. You can see that below then I specify the name of a new variable “crimr1” and then I tell the function I want that variable to equal (for each case) the division of the values in the variable “count” (number of crimes) by the variable “respop” (number of people residing in the area) and then we multiply the result of this division by 100,000 to obtain a rate expressed in those terms. Then we do likewise for the alternative measure of crime.

It should not be difficult for you to produce now a choropleth map like the one below. Clue: to change the colors for the fill of the polygons you can use the palette argument within the tm_fill calls. You can explore different palettes running the following code:

3.6.1 Homework 2

Reproduce the map below, you will need to include the code you used as part of your homework submission. Discuss the results. What are the most notable differences? Which denominator do you think is more appropriate (you will need to think about this quite carefully). Are you comparing like with like? Why? Why not? Could you make these comparisons more equivalent if you think you are not comparing like with like?*

## Some legend labels were too wide. These labels have been resized to 0.42, 0.42, 0.39. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
## Some legend labels were too wide. These labels have been resized to 0.42, 0.42, 0.42. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
## Legend labels were too wide. The labels have been resized to 0.42, 0.42, 0.42, 0.42, 0.42. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
## Some legend labels were too wide. These labels have been resized to 0.42, 0.42, 0.42, 0.39. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Once you have completed this activity, let’s explore your map with the crime rate using the usual residents as the denominator using the interactive way. Assuming you name that visualisation map5 you could use the following code.

## tmap mode set to interactive viewing
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.

You may find it useful to shift to the OpenStreetMap view by clicking in the box to the left, since it will give you a bit more contextual information than the default CartoDB basemap.

In the first lecture we spoke a bit about Open Street Map, but if you’re interested it’s definitely worth reading up on. As I mentioned, Open Street Map is a non-profit foundation whose aim is to support and enable the development of freely-reusable geospatial data, and relies heavily on volunteers participating in this project to map their local areas. You can have a look here for ongoing humanitarian projects, or read here about the mapping parties I told you about. At the very least though, in the spirit of open source and crowdsourcing, take a moment to appreciate that all these volunteers of people just like you have contributed to creating such a detailed geographical database of our world. That’s definitely something kinda cool to think about!

3.6.2 Homework 3

What areas of the city seem to have the largest concentration of crime? Why? Does deprivation help you to understand the patterns? What’s going on in the LSOA farthest to the South? Why does it look like a crime hotspot in this map and not in the one using the workday population?

3.7 More on small multiples and point pattern maps

Last week we showed you how to visualise point patterns using data from the Police UK website. One of the homeworks ask you to discuss the map you produced using leaflet in which you used type of crimes to colour your points. One of the problems with that map was that you had so many levels within that variable that it was very difficult to tell what was what. Even some of the colors in different categories were not that different from each other. That’s a typical situation where faceting or using small multiples would have been a better solution.

What we are going to do require our data to be stored as a spatial object -as it is the case with tmap. So first we need to turn our data into a spatial object. For this we will use the st_as_sf() function. The st_as_sf will return a sf object using the geographical coordinates we specify, below you can see we also specify the coordinate system we are using. Since we didn’t modify the tmap_mode from our last call we would still be running on the view rather than the plot format. We can go back to the plot format with a new tmap_mode call.

## tmap mode set to plotting

Here we are getting all the data points in the “crime_sf” object, which includes the whole of Greater Manchester. Also, since there are so many crimes, dots, it is hard to see patterns. We can add some transparency with the alpha argument as part of the tm_dots call.

Straight away you can start seeing some patterns, like the concentration of crimes in the town centres of all the different local authorities that conform Greater Manchester. This would be easier for you to see, if you are not familiar with the geography of Greater Manchester, if you place a basemap as a layer.

You can also use basemaps when working on the plot mode of tmap. We could for example use OpenStreet maps. We use the read_osm function from tmaptools package to do that, by passing “crime_sf” as an argument we will be bounding the scope of the basemap to the area covered by our point pattern of crimes. We will also need to load the OpenStreetMap package for this to work.

Remember that the first time you use a package you always have to install it, using the install.packages() function!

Now load the OpenStreetMap library into your working environment with the library() function:

Mac Users’s Note: So if you’re on a Mac, it’s possible that at this point you’ve experienced an error. It might look something like this:

Error: package or namespace load failed for ‘OpenStreetMap’: .onLoad failed in loadNamespace() for 'rJava', details: call: dyn.load(file, DLLpath = DLLpath, ...) error: unable to load shared object '/Library/Frameworks/R.framework/Versions/3.4/Resources/library/rJava/libs/rJava.so': dlopen(/Library/Frameworks/R.framework/Versions/3.4/Resources/library/rJava/libs/rJava.so, 6): Library not loaded: @rpath/libjvm.dylib Referenced from: /Library/Frameworks/R.framework/Versions/3.4/Resources/library/rJava/libs/rJava.so Reason: image not found

This is because of some drama between the rJava and your Mac OS’s dealing with Java. Do not despair though, the solution is simple. Open up your Terminal app. This is your command line for Mac. You might have used before for other things, or you might not. If you don’t know what I’m on about with Terminal have a look at this helpful tutorial.

Now, once you have Terminal open, all you have to do is copy and paste the code below, and press Enter to run it:

sudo ln -f -s $(/usr/libexec/java_home)/jre/lib/server/libjvm.dylib /usr/local/lib

After you press Enter, Terminal will ask you for your password. This is the password to your laptop. Type in your password, and hit enter again. Once that’s all done, you can go back to R, and you will have to load a package called rJava(). Like so:

Now you can again try loading the OpenStreetMap package.

Hopefully should all work smoothly now!

Back to everybody here!

You will also need another package called tmaptools. This gets installed when you installed tmap, but you still have to load it up.

We will then add the basemap layer using the qtm() function that is a function provided by tmap for quick plotting.

First, create an object, let’s call it gm_osm, by reading Open Street Map (OSM) data with the read_osm() function. This function reads and returns OSM tiles are read and returns as a spatial raster, or queries vectorized OSM data and returns as spatial polygons, lines, and/or points.

It might take a while to get this data, as it is quite a large set of data, being extracted from the open street map database. Just think how long it took all the volunteers to map these areas. Compared to that, the little waiting time to get this onto your computer for your map should be nothing!

Now you can use the qtm() function to draw a quick thematic map. Any ideas as to why it’s called qtm yet? Let’s try again. You can draw a Quick Thematic Map.

Again we use the tm_shape() and tm_dots() functions to specify the presentation of the map, with tm_shape() specifying the shape object (in this case crime_sf), and tm_dots() to draw symbols, including specifying the color, size, and shape of the symbols. In this case we adjust the transparency with the alpha parameter.

Let’s zoom in into Manchester city, for which we can use our manchester_lsoa map.

We can use the bb() function from tmaptols to create a bounding box around the University of Manchester, the width and height parameters I am using determine the degree of zoom in (I got these experimenting with different ones until I got the right zooming around the University location). Once you have this bounding box you can pass it as an argument to the read_osm() function that will look for the basemap around that location.

Now that we have our basemap we can run our small multiples. Because we have a variable that determines the types we can use a different way to tmap_arrange() explained above, we can essentially ask tmap to create a map for each of the levels in our organising variable (in this case Crime.type). So instead of tmap_arrange() that requires the creation of each map, when each map simply represents different levels of an organising variable we can simplify the syntax using tm_facets() and within this function we specify as the first argument the name of the variable that has the different categories we want to map out. The second argument you see below free.coords set to FALSE simply ensures that the map gets bounded to the basemap, if you want to see what happens if you change it, just set it to TRUE instead.

We could do some further tweaking around for ensuring things look a bit neater. But we have covered a lot of ground today, and you should all give yourself a congratulatory “well done”.